Sensitivity analysis: Functional Entity grouping traits
Compiled on Wed Oct 15 14:38:54 2025 by ohara
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE,
fig.height = 4, fig.width = 7)
library(terra)
library(sf)
library(oharac)
library(data.table)
library(tidyverse)
library(here)
library(cowplot)
source(here('common_fxns.R'))1 Summary
Sensitivity analysis to determine how drastically scores change with different values of functional traits.
In script 4a, we calculated sensitivity around functional diversity metrics and per-stressor impacts. In 4b, we aggregate results across all stressors for each randomized iteration, to determine the distribution of mean values for cumulative impacts. Here, we map them out and make cool figures.
2 Methods
For each randomized trait, plot the variation in functional diversity metrics for the sampled points, with a nearest-neighbor interpolation to help spot regional variation. We will examine several stats: absolute and relative difference (diff / mean), 95% confidence interval, coefficient of variance (std dev divided by mean)
ocean_r <- rast(here('_spatial/ocean_area_mol.tif'))
cell_xy <- ocean_r %>%
setValues(1:ncell(.)) %>%
setNames('cell_id') %>%
as.data.frame(xy = TRUE)fd_fs <- list.files(here_anx('sens_analysis/funct_div_summary'),
pattern = '.csv', full.names = TRUE)
fd_ts <- basename(fd_fs) %>% str_remove_all('.+_shuffle_|.csv')
fd_df <- lapply(fd_fs, fread) %>%
setNames(fd_ts) %>%
rbindlist(idcol = 't') %>%
mutate(n_fe_diff = n_fe_mean - n_fe,
fv_diff = fv_mean - fv,
n_fe_ci = (n_fe_95hi - n_fe_95lo),
fv_ci = (fv_95hi - fv_95lo),
n_fe_cv = n_fe_sdev / n_fe_mean,
fv_cv = fv_sdev / fv_mean) %>%
rename(n_fe_value = n_fe, fv_value = fv) %>%
pivot_longer(cols = -c(cell_id, t, nspp),
names_to = 'var', values_to = 'val') %>%
mutate(param = str_extract(var, 'n_fe|fv'),
stat = str_extract(var, 'mean|sdev|ci|cv|diff|95..|value')) %>%
select(-var) %>%
left_join(cell_xy, by = 'cell_id')2.1 Functional diversity metrics
2.1.1 Maps of functional vulnerability
How does functional vulnerability vary as traits are randomized?
land_sf <- read_sf(here('_spatial/ne_10m_land/ne_10m_land_no_casp.shp')) %>%
st_transform(st_crs(ocean_r))
globe_bbox <- rbind(c(-180, -90), c(-180, 90),
c(180, 90), c(180, -90), c(-180, -90))
globe_border <- st_polygon(list(globe_bbox)) %>%
st_sfc(crs = 4326) %>%
st_sf(data.frame(rgn = 'globe', geom = .)) %>%
smoothr::densify(max_distance = 0.5) %>%
st_transform(crs = crs(ocean_r))map_param <- function(df, s, p, r = ocean_r, sig = FALSE, thresh = NULL) {
val_mtx <- df %>%
filter(stat == s & param == p) %>%
select(x, y, z = val) %>%
as.matrix()
val_rast <- interpNear(r, val_mtx, radius = 500000) %>%
mask(r)
if(sig) {
sig_mtx <- df %>%
filter(param == p & stat %in% c('95hi', '95lo', 'value')) %>%
pivot_wider(names_from = stat, values_from = val) %>%
mutate(sig = ifelse(between(value, `95lo`, `95hi`), 0, 1)) %>%
select(x, y, z = sig) %>%
as.matrix()
sig_rast <- interpNear(r, sig_mtx, radius = 50000) %>%
mask(r)
val_rast[sig_rast == 0] <- NA
}
if(!is.null(thresh)) val_rast[abs(val_rast) < thresh] <- NA
return(val_rast)
}
finish_map <- function(p, div = TRUE) {
### div is whether the palette should be a diverging palette (TRUE) or
### sequential (FALSE)
p <- p +
### continents:
geom_sf(data = land_sf,
fill = 'grey96', color = 'grey40',
size = .10) +
geom_sf(data = globe_border,
fill = NA, color = 'grey70',
size = .1) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
facet_wrap(~ trait) +
theme_void()
if(div) {
### diverging palette
p <- p +
scale_fill_gradient2(low = '#d01c8b', mid = '#ffffdf', high = '#4dac26')
} else {
### sequential palette
p <- p + scale_fill_viridis_c()
}
return(p)
}fv_fig_f <- '_figs/fv_diff_by_trait_map.png'
if(!file.exists(fv_fig_f)) {
fv_map_list <- vector('list', length = length(fd_ts)) %>%
setNames(fd_ts)
for(trait in fd_ts) {
# trait <- fd_ts[1]
fv_map_list[[trait]] <- map_param(df = fd_df %>% filter(t == trait),
p = 'fv', s = 'diff')
}
fv_map_stack <- rast(fv_map_list)
mean_sd_vals <- as.data.frame(fv_map_stack) %>%
summarize(across(everything(), .fns = list(mean = mean, sd = sd), .names = '{.col}_{.fn}'))
# len_mean len_sd mob_mean mob_sd trp_mean trp_sd wcol_mean wcol_sd
# 1 0.03633012 0.04194202 0.03743112 0.04292291 0.05350717 0.04770103 0.02206068 0.02724206
fv_map_df <- as.data.frame(fv_map_stack, xy = TRUE) %>%
pivot_longer(-c(x, y), names_to = 'trait', values_to = 'delta_fv') %>%
mutate(trait = case_when(trait == 'len' ~ 'body length',
trait == 'mob' ~ 'adult mobility',
trait == 'trp' ~ 'trophic level',
trait == 'wcol' ~ 'water column position'))
p <- ggplot() +
geom_raster(data = fv_map_df, aes(x, y, fill = delta_fv)) +
labs(fill = 'ΔFV')
p <- finish_map(p)
ggsave(fv_fig_f, width = 6.5, height = 3, dpi = 300)
}
knitr::include_graphics(fv_fig_f)
fv2_fig_f <- '_figs/fv_cv_by_trait_map.png'
if(!file.exists(fv2_fig_f)) {
fv_map_list <- vector('list', length = length(fd_ts)) %>%
setNames(fd_ts)
for(trait in fd_ts) {
fv_map_list[[trait]] <- map_param(df = fd_df %>% filter(t == trait),
p = 'fv', s = 'cv')
}
fv_map_stack <- rast(fv_map_list)
fv_map_df <- as.data.frame(fv_map_stack, xy = TRUE) %>%
pivot_longer(-c(x, y), names_to = 'trait', values_to = 'cv_fv') %>%
mutate(trait = case_when(trait == 'len' ~ 'body length',
trait == 'mob' ~ 'adult mobility',
trait == 'trp' ~ 'trophic level',
trait == 'wcol' ~ 'water column position'))
p <- ggplot() +
geom_raster(data = fv_map_df, aes(x, y, fill = cv_fv)) +
labs(fill = 'CV(FV)')
p <- finish_map(p, div = FALSE)
ggsave(fv2_fig_f, width = 6.5, height = 3, dpi = 300)
}
knitr::include_graphics(fv2_fig_f)
2.1.2 Plots of functional diversity parameters vs. number of species
How do functional diversity metrics (FV and number of FEs) vary as traits are randomized, relative to the number of species in a given pixel?
plot_f <- '_figs/fv_stat_by_trait_scatter.png'
if(!file.exists(plot_f)) {
fv_plot_df <- fd_df %>%
filter(stat %in% c('mean', 'cv', 'diff')) %>%
filter(param == 'fv') %>%
mutate(metric = paste(stat, param)) %>%
mutate(metric = str_replace(metric, 'fv', '(FV)'),
metric = str_replace(metric, 'mean ', 'mean'),
metric = str_replace(metric, 'cv ', 'CV'),
metric = str_replace(metric, 'diff ', 'Δ')) %>%
mutate(t = case_when(t == 'len' ~ 'body length',
t == 'mob' ~ 'adult mobility',
t == 'trp' ~ 'trophic level',
t == 'wcol' ~ 'water column position')) %>%
sample_frac(.1)
p <- ggplot(fv_plot_df, aes(x = nspp, y = val)) +
geom_point(alpha = .2) +
theme_ohara() +
facet_grid(metric ~ t, scales = 'free_y') +
labs(x = 'Species richness') +
theme(axis.title.y = element_blank())
ggsave(plot_f, width = 6, height = 3, dpi = 300)
}
knitr::include_graphics(plot_f)
plot_f <- '_figs/n_fe_stat_by_trait_scatter.png'
if(!file.exists(plot_f)) {
fe_plot_df <- fd_df %>%
filter(stat %in% c('mean', 'cv', 'diff')) %>%
filter(param == 'n_fe') %>%
mutate(metric = paste(stat, param)) %>%
mutate(metric = str_replace(metric, 'n_fe', '(# FE)'),
metric = str_replace(metric, 'mean ', 'mean'),
metric = str_replace(metric, 'cv ', 'CV'),
metric = str_replace(metric, 'diff ', 'Δ')) %>%
mutate(t = case_when(t == 'len' ~ 'body length',
t == 'mob' ~ 'adult mobility',
t == 'trp' ~ 'trophic level',
t == 'wcol' ~ 'water column position')) %>%
sample_frac(.1)
p <- ggplot(fe_plot_df, aes(x = nspp, y = val)) +
geom_point(alpha = .2) +
theme_ohara() +
facet_grid(metric ~ t, scales = 'free_y') +
labs(x = 'Species richness') +
theme(axis.title.y = element_blank())
ggsave(plot_f, width = 6, height = 3, dpi = 300)
}
knitr::include_graphics(plot_f)
2.2 Cumulative impact metrics
2.2.1 Maps of cumulative impact variation
How does cumulative impact vary as traits are randomized?
chi_fs <- list.files(here_anx('sens_analysis/impact_summary'),
pattern = 'chi_.+.csv', full.names = TRUE)
chi_ts <- basename(chi_fs) %>% str_remove_all('.+_shuffle_|.csv')
chi_true_df <- rast(here('_output/cumulative_impact_maps/chi_funct_entity.tif')) %>%
as.data.frame(xy = TRUE) %>%
oharac::dt_join(cell_xy, by = c('x', 'y'), type = 'left') %>%
filter(cell_id %in% fd_df$cell_id)
nspp_df <- fd_df %>%
select(cell_id, nspp) %>%
distinct()
chi_df <- lapply(chi_fs, fread) %>%
setNames(chi_ts) %>%
rbindlist(idcol = 't') %>%
left_join(chi_true_df, by = 'cell_id') %>%
left_join(nspp_df, by = 'cell_id') %>%
mutate(chi_diff = chi_mean - chi_fe,
chi_ci = (chi_95hi - chi_95lo),
chi_cv = chi_sdev / chi_mean) %>%
select(cell_id, t, nspp,
chi_mean, chi_diff,
chi_ci, chi_cv) %>%
pivot_longer(cols = -c(cell_id, t, nspp),
names_to = 'var', values_to = 'val') %>%
mutate(param = 'chi',
stat = str_extract(var, 'mean|sdev|ci|cv|diff')) %>%
select(-var) %>%
left_join(cell_xy, by = 'cell_id')chi_fig_f <- '_figs/chi_diff_by_trait_map.png'
if(!file.exists(chi_fig_f)) {
chi_map_list <- vector('list', length = length(chi_ts)) %>%
setNames(chi_ts)
for(trait in fd_ts) {
chi_map_list[[trait]] <- map_param(df = chi_df %>% filter(t == trait),
p = 'chi', s = 'diff')
}
chi_map_stack <- rast(chi_map_list)
mean_sd_vals <- as.data.frame(chi_map_stack) %>%
summarize(across(everything(), .fns = list(mean = mean, sd = sd), .names = '{.col}_{.fn}'))
# len_mean len_sd mob_mean mob_sd trp_mean trp_sd wcol_mean wcol_sd
# 1 -0.00801914 0.01990471 0.0008610911 0.01656777 0.003189031 0.01596556 -0.006297416 0.02280766
chi_map_df <- as.data.frame(chi_map_stack, xy = TRUE) %>%
pivot_longer(-c(x, y), names_to = 'trait', values_to = 'delta_chi') %>%
mutate(trait = case_when(trait == 'len' ~ 'body length',
trait == 'mob' ~ 'adult mobility',
trait == 'trp' ~ 'trophic level',
trait == 'wcol' ~ 'water column position'))
p <- ggplot() +
geom_raster(data = chi_map_df, aes(x, y, fill = delta_chi)) +
labs(fill = 'ΔCHI')
p <- finish_map(p)
ggsave(chi_fig_f, width = 6.5, height = 3, dpi = 300)
}
knitr::include_graphics(chi_fig_f)
2.2.2 Plots of CHI stats vs. number of species
How do CHI stats vary as traits are randomized, relative to the number of species in a given pixel?
plot_f <- '_figs/chi_stat_by_trait_scatter.png'
if(!file.exists(plot_f)) {
chi_plot_df <- chi_df %>%
filter(stat %in% c('mean', 'cv', 'diff')) %>%
mutate(metric = paste(stat, param)) %>%
mutate(metric = str_replace(metric, 'chi', '(CHI)'),
metric = str_replace(metric, 'mean ', 'mean'),
metric = str_replace(metric, 'cv ', 'CV'),
metric = str_replace(metric, 'diff ', 'Δ')) %>%
mutate(t = case_when(t == 'len' ~ 'body length',
t == 'mob' ~ 'adult mobility',
t == 'trp' ~ 'trophic level',
t == 'wcol' ~ 'water column position')) %>%
sample_frac(.1)
p <- ggplot(chi_plot_df, aes(x = nspp, y = val)) +
geom_point(alpha = .2) +
theme_ohara() +
facet_grid(metric ~ t, scales = 'free_y') +
labs(x = 'Species richness') +
theme(axis.title.y = element_blank())
ggsave(plot_f, width = 6, height = 3, dpi = 300)
}
knitr::include_graphics(plot_f)
3 Fig for SI
Combine the FV variation and CHI variation into a large panel.
p <- ggdraw() +
draw_image(fv_fig_f, x = 0.0, y = 0.51, width = 1.0, height = 0.49) +
draw_image(chi_fig_f, x = 0.0, y = 0.0, width = 1.0, height = 0.49) +
draw_label('A.', x = .07, y = .995, hjust = 0, vjust = 1,
size = 7, color = 'black', fontfamily = 'Helvetica', fontface = 'bold') +
draw_label('B.', x = 0.5, y = .995, hjust = 0, vjust = 1,
size = 7, color = 'black', fontfamily = 'Helvetica', fontface = 'bold') +
draw_label('C.', x = .07, y = .748, hjust = 0, vjust = 1,
size = 7, color = 'black', fontfamily = 'Helvetica', fontface = 'bold') +
draw_label('D.', x = 0.5, y = .748, hjust = 0, vjust = 1,
size = 7, color = 'black', fontfamily = 'Helvetica', fontface = 'bold') +
draw_label('E.', x = .07, y = .485, hjust = 0, vjust = 1,
size = 7, color = 'black', fontfamily = 'Helvetica', fontface = 'bold') +
draw_label('F.', x = 0.5, y = .485, hjust = 0, vjust = 1,
size = 7, color = 'black', fontfamily = 'Helvetica', fontface = 'bold') +
draw_label('G.', x = .07, y = .235, hjust = 0, vjust = 1,
size = 7, color = 'black', fontfamily = 'Helvetica', fontface = 'bold') +
draw_label('H.', x = 0.5, y = .235, hjust = 0, vjust = 1,
size = 7, color = 'black', fontfamily = 'Helvetica', fontface = 'bold')
si_fig_file <- here('5_ms_figs/figS6_sens_analysis.png')
ggsave(si_fig_file, height = 11.5, width = 12, units = 'cm', dpi = 300)
knitr::include_graphics(si_fig_file)